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Abstract 

The global energy fluctuations of a low density gas granular gas in the homogeneous cooling state 
near its clustering instability are studied by means of molecular dynamics simulations. The relative 
dispersion of the fluctuations is shown to exhibit a power-law divergent behavior. Moreover, the 
probability distribution of the fluctuations presents data collapse as the system approaches the 
instability, for different values of the inelasticity. The function describing the collapse turns out to 
be the same as the one found in several molecular equilibrium and non-equilibrium systems, except 
for the change in the sign of the fluctuations. 

PACS numbers: PACS Numbers: 45.70.-n,51.10.+y,05.20.Dd 



A granular system is an assembly of macroscopic particles dissipating their energy through 
inelastic collisions. In the rapid granular flow regime, the grains move freely and indepen- 
dently between collisions and the system is often referred to as a granular gas. In this regime, 
the behavior of a granular system resembles that of a molecular fluid, although with many 
significant differences . These discrepancies have been extensively illustrated by using the- 
oretical, experimental, and also particle simulation methods In this paper, evidence will 
be reported of a rather surprising analogy between granular gases and molecular systems. 
It will be shown that the statistical behavior of a freely evolving granular system near the 
clustering instability has strong similarities with that of several molecular model systems, 
including equilibrium systems in the critical region. More precisely, the second moment of 
the relative fluctuations of the global energy diverges with a given critical exponent. Besides, 
and even more remarkably, the energy fluctuations obey a scaling law that, when properly 
expressed, is the same as the one found in some correlated equilibrium and non-equilibrium 
molecular systems 0, 0, 0] . 

A widely used simple model for granular gases is a system of smooth inelastic hard 
spheres {d = 3) or disks {d = 2) of mass m and diameter a. The inelasticity of collisions 
is characterized by a constant coefficient of normal restitution a. As a consequence of the 
energy dissipation in collisions, isolated granular gases do not exhibit any homogeneous time- 
independent state. The simplest state is the so-called homogeneous cooling state (HCS) 
with vanishing flow field and a monotonically decreasing temperature Thcs{t), obeying the 
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equation dtTtcsit) = -Chcs{Thcs)Thcs{t), where Chcs{Thcs) oc T^^ is the cooling rate. At a 
microscopic level, it has been postulated that the time dependence of the ensemble describing 
this state occurs only through the scaling of the velocities with the thermal velocity, that 
is proportional to T;[£^(t), and the corresponding normalization 1?^. The HCS is unstable 
against long wavelength spatial perturbations (sj, leading to the spontaneous formation of 
velocity vortices and density clusters. This clustering instability is accurately predicted by 



a linear stability analysis_p: 
transversal shear mode 



the hydrodynamic equations, that shows that it is driven by the 



y, ^. A critical length Lc is identified, so that the system becomes 
unstable when its linear size L exceeds L^. Its value has been determined for a low density 
granular gas described by the (inelastic) Boltzmann equation, and it is given by Q| 
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where n is the number of particles density, C*('^) = VoChcs/nkBThcs, and ri*{a) = 
viThcs) / VoiThcs) , with kB being the Boltzmann constant, rj the shear viscosity of the granular 
gas, and rjQ its elastic limit. The explicit expressions of (* and rj* are given in In the 
elastic limit a = 1, vanishes and rj* = l.The accuracy of this prediction has been verified 



ll|. 



by direct Monte Carlo simulation of the Boltzmann equation 

Another important consequence of the inelasticity of collisions, largely unexplored, is 
;he presence of an intrinsic noise in the macroscopic description of the system. Recently 
1^ . it has been shown that the global energy of a dilute granular gas in the HCS exhibits 
fluctuations and time-correlation properties which are caused by the energy dissipation in 
collisions and, at a mesoscopic level, by the presence of velocity correlations. Although the 
correlations increase as the inelasticity increases, they remain relatively small over all the 
range of values of a studied {a > 0.6). Here it will be shown that the noise is amplified, 
and in fact it diverges, at the threshold of the clustering instability. 

We have performed molecular dynamics (MD) simulations of a freely evolving system of 
N inelastic hard disks (rf = 2), using the steady-state method, that is based on an exact 
mapping of the HCS onto a steady state. This property follows from a change in the time 
scale, as discussed in detail in [l3[. The steady representation of the HCS removes the 
difficulties associated with the rapid cooling of the fluid leading to numerical inaccuracies 
very soon. With this method, the trajectories can be followed for an arbitrary time. The 
system considered has been a square box of size L with periodic boundary conditions. The 
density in all the results reported in the following is n = 0.02(T~^, so the values of L and 
N were consistently changed in the different simulations. The above value of the density 
has been found to be low enough as to guarantee that the average behavior of the system is 
accurately described by the Boltzmann equation, even in the unstable region of the granular 
gas ^|. Two different values of the restitution coefficient have been considered, a = 0.8 and 
a = 0.9, for which Eq. ^ gives Lc — 304(T and Lc — 4:13a, respectively. We have studied 
the behavior of the global properties of the system when the critical size is approached 
from below by increasing the value of L at constant a (and n). In all the simulations, we 
have checked that the system was actually in the HCS by monitoring the local velocity and 
density fluctuations. 

The first quantity we have investigated is the reduced cooling rate (* as a function of 
the linear size of the system, at constant a. The simulations show that it monotonically 
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FIG. 1: Dependence of the scaled cooling rate C*(-^^) on the "distance" to the instability point 
6L = I — L/Lc for a system of inelastic hard disks with a = 0.9. The solid line is a fit to the 
power-law behavior discussed in the text. 

decreases very slowly as L approaches its critical value. The decay is imperceptible on the 

scale used in Fig. Q where ln(*~'^ is plotted as a function of 6L = (Lc — L)/Lc for a = 0.9. 

Here and in the following, the theoretical prediction for Lc given by Eq. has been used. 

To make the above dependence on L explicit, we have also plotted ln[(*~^{L) — C*^^(0)], 

where C*(0) is the constant asymptotic value of the reduced cooling rate far below the 

transition point, also obtained from the simulations. The results indicate a behavior of the 

form (*~^{L) ~ (*~^{0) + A(^6L \ with ~ 0.02. Similar results were found for a = 0.8. If 

this behavior persists until the instability point, it would imply that the cooling rate tends 

to vanish as (^L/A^)^/^. Nevertheless, since in the observed region C*{L) is dominated by 

its asymptotic part C*(0), it is difficult to make any definite statement on this point. 

Next, let us consider the second moment of the fiuctuations of the total energy E, 

, _ {E\t)) - {E{t)f 
{Eit))^ ' 

where the angular brackets denote average over the ensemble generated by a number of 
trajectories [l^. The above quantity does not depend on time in the HCS, due to the scaling 
property of its distribution function. The results, again for a = 0.9, are given in Fig. |21 
For L <^ Lc, Na%{L) is practically constant, independent of L, indicating that aE oc A^"^/^, 
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FIG. 2: Scaled second moment of the total energy fluctuations (7 ^ SiS 3j function of the "distance" 
6L to the instability point for a = 0.9. Note that the second moment has been multiplied by the 
number of particles A^. The solid line is the fit to the power law defined in the text. 

as expected. In fact, its value in this region, that we will denote by Na'^{0), is accurately 



predicted by the result derived in ref. 



12| by using kinetic theory methods. Although it is a 



function of a, increasing as a decreases, it remains rather small, at least for not too inelastic 
systems. On the other hand, when L approaches its critical value, the energy fluctuations 
grow very fast. This is clearer seen by studying the quantity [cM-^) " '^e{^)]^ which is 
also plotted in the same figure. It follows that close enough to the instability point, the 
second moment of the total energy fluctuations is accurately described by the critical law 
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where Ae — Q x 10~^ is a critical amplitude. Upon writing the above equation, we have 
used that near the instability point we can approximate N ^ Nc = nL^. The same behavior 
is found for a = 0.8, but with the amplitude Ae — 1.5 x 10"'^. Therefore, the simulation 
results indicate a divergent critical behavior of the fluctuations, with a critical exponent —1 
and an amplitude that depends on the value of the coefficient of restitution a. Of course, 
the validity of the above relies on the assumption that the observed behavior remains the 
same until the system is asymptotically close to the instability. 

To analyze in more detail the structure of the energy fluctuations, we make the unjustified 



assumption {a priori) that the probabihty distribution function (PDF) Pi{5E) of the relative 
energy fluctuations 5E = {E(t) — {E{t)))/ {E(t)) in the vicinity of the clustering instability 
verifies the scaling relation 



where / is a scaling function. It is moreover assumed that all the dependence of / on a and 
L occurs through aE- Consider then the function 



All the L and a dependence has been eliminated on the right hand side of the above expres- 
sion when considered as a function of 6E/aE- Thus as long as the assumed scaling holds 
near the clustering instability, the data for systems with different sizes should tend to col- 
lapse onto each other as L approaches Lc. Moreover, the collapse must occur on a function 
that does not depend on a. To check this strong prediction, we have computed the PDF 
of the fluctuations in the MD simulations. For each trajectory of a system, the trace of SE 
has been partitioned into nonoverlapping bins and the frequency distribution has been built 
up using all the trajectories corresponding to the same values of L and a. The resulting 
normalized distributions are plotted in Figs. ElandlU for a = 0.8 and a = 0.9, respectively. 
In each case, results for different system sizes are given. Notice that, for reasons that will be 
explained below, what has been actually plotted is aEPii'SE). The data are seen to fall 
quite closely on top of each other, as implied by the scaling law (0), over more than three 
orders of magnitude, specially when attention is restricted to the results for the two largest 
systems. Moreover, the fluctuations are highly non-Gaussian and asymmetric around the 
mean value. This is a manifestation of the presence of a correlation length of the order of 
the size of the system, so that it can not be divided into statistical independent mesoscopic 
regions and, therefore, there is no reason to expect the fluctuations of global quantities to 
be Gaussian. 

The shape of PDF's in Figs. El and ID looks quite similar to the functional form found 
in other equilibrium and non-equilibrium systems 0, 0]. To make this statement more 
precise, we have considered the approximated form of the PDF for the fluctuations of the 
global magnetization in the two-dimensional XY model in the harmonic approximation, for 
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FIG. 3: Probability density function of the relative total energy fluctuations aEPii—^E) for a 
system of inelastic hard disks with a = 0.8. The symbols are from MD simulations and the solid 
line is Eq. ®. 



temperatures well below the Kosterlitz-Thouless transition 

Ii{y) = K (e^-^y , x = biy~s), a = 7r/2. (6) 

The high y region of this function is approximately Gaussian, while it presents an exponential 
tail for large negative values. A quite large number of non-equilibrium systems exhibiting 
self-organized criticality has been found to show data collapse with a PDF very similar 
to n(y) Quite remarkably, the same kind of behavior had been previously observed 
in experiments with confined turbulent flows 2. In the present context, in principle it 
should be Il{y) = crEPiiy) with y = 6E, but in order to get a good agrement with the 
MD results, we have to change 6E by —6E in the identification of y. The values of the 
three parameters in Eq. © follow from the normalization, the zero mean, and unit variance 
conditions, with the results K = 2.14, b = 0.938, and s = 0.374. Therefore, Eq. © has 
no fitting parameters. The function Il{y) has been also plotted in Figs. El and EJ and a 
surprisingly good agreement with the simulation data is found. The change in the sign of 
the fluctuations, i.e. the fact that the energy fluctuations in a granular gas are described by 
a function that is the symmetric with respect to the origin of that for ordinary systems, may 
be due to the dissipative character of granular systems. In the HCS the gas is continuously 








FIG. 4: The same as in Fig. |21 but for a = 0.9. 

dissipating energy, while in molecular systems energy must be continuously supplied in 
order to keep them in a nonequilibrium steady state. Let us mention that large fluctuations 
around the threshold of a symmetry breaking instability have also been observed recently 
in granular systems driven by a thermal sidewall at zero gravity To our knowledge, 

whether these fluctuations can be scaled in a similar way to the one described here has not 
been investigated yet. 

The natural question prompted for the results reported in this work is whether the agree- 
ment of the fluctuation spectra of such a variety of systems is telling us something about the 
intrinsic behavior of a quite general class of systems. If this is the case, how this class can be 
characterized a priori, and the divergent behavior (jH)) and the specific scaling form (jHl), or 
another one very close to it, can be derived on the ground of general arguments. Specifically 
in the context of granular systems, how can the above results be obtained starting from a 
microscopic description of the system? Are they captured by a fluctuating hydrodynamics 
theory? Is there an underlying hyperscaling in the sense discussed in jsl? Finally, it must 
be stressed that we have restricted ourselves here to two-dimensional low-density granular 
gases. Although it seems that the same kind of results can be expected at higher densities 
and also for three-dimensional systems (perhaps with a different critical exponent), this is 
something to be verified. These points are presently under study and some results will be 
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reported elsewhere. 
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